<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>dlarrv.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="DLARRV.1"></a><a href="dlarrv.f.html#DLARRV.1">DLARRV</a>( N, VL, VU, D, L, PIVMIN,
     $                   ISPLIT, M, DOL, DOU, MINRGP,
     $                   RTOL1, RTOL2, W, WERR, WGAP,
     $                   IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
     $                   WORK, IWORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK auxiliary routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      INTEGER            DOL, DOU, INFO, LDZ, M, N
      DOUBLE PRECISION   MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      INTEGER            IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
     $                   ISUPPZ( * ), IWORK( * )
      DOUBLE PRECISION   D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
     $                   WGAP( * ), WORK( * )
      DOUBLE PRECISION  Z( LDZ, * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="DLARRV.26"></a><a href="dlarrv.f.html#DLARRV.1">DLARRV</a> computes the eigenvectors of the tridiagonal matrix
</span><span class="comment">*</span><span class="comment">  T = L D L^T given L, D and APPROXIMATIONS to the eigenvalues of L D L^T.
</span><span class="comment">*</span><span class="comment">  The input eigenvalues should have been computed by <a name="DLARRE.28"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The order of the matrix.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  VL      (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">  VU      (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">          Lower and upper bounds of the interval that contains the desired
</span><span class="comment">*</span><span class="comment">          eigenvalues. VL &lt; VU. Needed to compute gaps on the left or right
</span><span class="comment">*</span><span class="comment">          end of the extremal eigenvalues in the desired RANGE.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D       (input/output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, the N diagonal elements of the diagonal matrix D.
</span><span class="comment">*</span><span class="comment">          On exit, D may be overwritten.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  L       (input/output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, the (N-1) subdiagonal elements of the unit
</span><span class="comment">*</span><span class="comment">          bidiagonal matrix L are in elements 1 to N-1 of L
</span><span class="comment">*</span><span class="comment">          (if the matrix is not splitted.) At the end of each block
</span><span class="comment">*</span><span class="comment">          is stored the corresponding shift as given by <a name="DLARRE.50"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a>.
</span><span class="comment">*</span><span class="comment">          On exit, L is overwritten.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  PIVMIN  (in) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">          The minimum pivot allowed in the Sturm sequence.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  ISPLIT  (input) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The splitting points, at which T breaks up into blocks.
</span><span class="comment">*</span><span class="comment">          The first block consists of rows/columns 1 to
</span><span class="comment">*</span><span class="comment">          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
</span><span class="comment">*</span><span class="comment">          through ISPLIT( 2 ), etc.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  M       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The total number of input eigenvalues.  0 &lt;= M &lt;= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  DOL     (input) INTEGER
</span><span class="comment">*</span><span class="comment">  DOU     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          If the user wants to compute only selected eigenvectors from all
</span><span class="comment">*</span><span class="comment">          the eigenvalues supplied, he can specify an index range DOL:DOU.
</span><span class="comment">*</span><span class="comment">          Or else the setting DOL=1, DOU=M should be applied.
</span><span class="comment">*</span><span class="comment">          Note that DOL and DOU refer to the order in which the eigenvalues
</span><span class="comment">*</span><span class="comment">          are stored in W.
</span><span class="comment">*</span><span class="comment">          If the user wants to compute only selected eigenpairs, then
</span><span class="comment">*</span><span class="comment">          the columns DOL-1 to DOU+1 of the eigenvector space Z contain the
</span><span class="comment">*</span><span class="comment">          computed eigenvectors. All other columns of Z are set to zero.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  MINRGP  (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RTOL1   (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">  RTOL2   (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">           Parameters for bisection.
</span><span class="comment">*</span><span class="comment">           An interval [LEFT,RIGHT] has converged if
</span><span class="comment">*</span><span class="comment">           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  W       (input/output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The first M elements of W contain the APPROXIMATE eigenvalues for
</span><span class="comment">*</span><span class="comment">          which eigenvectors are to be computed.  The eigenvalues
</span><span class="comment">*</span><span class="comment">          should be grouped by split-off block and ordered from
</span><span class="comment">*</span><span class="comment">          smallest to largest within the block ( The output array
</span><span class="comment">*</span><span class="comment">          W from <a name="DLARRE.89"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a> is expected here ). Furthermore, they are with
</span><span class="comment">*</span><span class="comment">          respect to the shift of the corresponding root representation
</span><span class="comment">*</span><span class="comment">          for their block. On exit, W holds the eigenvalues of the
</span><span class="comment">*</span><span class="comment">          UNshifted matrix.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WERR    (input/output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The first M elements contain the semiwidth of the uncertainty
</span><span class="comment">*</span><span class="comment">          interval of the corresponding eigenvalue in W
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WGAP    (input/output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The separation from the right neighbor eigenvalue in W.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IBLOCK  (input) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The indices of the blocks (submatrices) associated with the
</span><span class="comment">*</span><span class="comment">          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
</span><span class="comment">*</span><span class="comment">          W(i) belongs to the first block from the top, =2 if W(i)
</span><span class="comment">*</span><span class="comment">          belongs to the second block, etc.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INDEXW  (input) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The indices of the eigenvalues within each block (submatrix);
</span><span class="comment">*</span><span class="comment">          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
</span><span class="comment">*</span><span class="comment">          i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  GERS    (input) DOUBLE PRECISION array, dimension (2*N)
</span><span class="comment">*</span><span class="comment">          The N Gerschgorin intervals (the i-th Gerschgorin interval
</span><span class="comment">*</span><span class="comment">          is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should
</span><span class="comment">*</span><span class="comment">          be computed from the original UNshifted matrix.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Z       (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
</span><span class="comment">*</span><span class="comment">          If INFO = 0, the first M columns of Z contain the
</span><span class="comment">*</span><span class="comment">          orthonormal eigenvectors of the matrix T
</span><span class="comment">*</span><span class="comment">          corresponding to the input eigenvalues, with the i-th
</span><span class="comment">*</span><span class="comment">          column of Z holding the eigenvector associated with W(i).
</span><span class="comment">*</span><span class="comment">          Note: the user must ensure that at least max(1,M) columns are
</span><span class="comment">*</span><span class="comment">          supplied in the array Z.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDZ     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array Z.  LDZ &gt;= 1, and if
</span><span class="comment">*</span><span class="comment">          JOBZ = 'V', LDZ &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  ISUPPZ  (output) INTEGER array, dimension ( 2*max(1,M) )
</span><span class="comment">*</span><span class="comment">          The support of the eigenvectors in Z, i.e., the indices
</span><span class="comment">*</span><span class="comment">          indicating the nonzero elements in Z. The I-th eigenvector
</span><span class="comment">*</span><span class="comment">          is nonzero only in elements ISUPPZ( 2*I-1 ) through
</span><span class="comment">*</span><span class="comment">          ISUPPZ( 2*I ).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace) DOUBLE PRECISION array, dimension (12*N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IWORK   (workspace) INTEGER array, dimension (7*N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          &gt; 0:  A problem occured in <a name="DLARRV.142"></a><a href="dlarrv.f.html#DLARRV.1">DLARRV</a>.
</span><span class="comment">*</span><span class="comment">          &lt; 0:  One of the called subroutines signaled an internal problem.
</span><span class="comment">*</span><span class="comment">                Needs inspection of the corresponding parameter IINFO
</span><span class="comment">*</span><span class="comment">                for further information.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          =-1:  Problem in <a name="DLARRB.147"></a><a href="dlarrb.f.html#DLARRB.1">DLARRB</a> when refining a child's eigenvalues.
</span><span class="comment">*</span><span class="comment">          =-2:  Problem in <a name="DLARRF.148"></a><a href="dlarrf.f.html#DLARRF.1">DLARRF</a> when computing the RRR of a child.
</span><span class="comment">*</span><span class="comment">                When a child is inside a tight cluster, it can be difficult
</span><span class="comment">*</span><span class="comment">                to find an RRR. A partial remedy from the user's point of
</span><span class="comment">*</span><span class="comment">                view is to make the parameter MINRGP smaller and recompile.
</span><span class="comment">*</span><span class="comment">                However, as the orthogonality of the computed vectors is
</span><span class="comment">*</span><span class="comment">                proportional to 1/MINRGP, the user should be aware that
</span><span class="comment">*</span><span class="comment">                he might be trading in precision when he decreases MINRGP.
</span><span class="comment">*</span><span class="comment">          =-3:  Problem in <a name="DLARRB.155"></a><a href="dlarrb.f.html#DLARRB.1">DLARRB</a> when refining a single eigenvalue
</span><span class="comment">*</span><span class="comment">                after the Rayleigh correction was rejected.
</span><span class="comment">*</span><span class="comment">          = 5:  The Rayleigh Quotient Iteration failed to converge to
</span><span class="comment">*</span><span class="comment">                full accuracy in MAXITR steps.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Based on contributions by
</span><span class="comment">*</span><span class="comment">     Beresford Parlett, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment">     Jim Demmel, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment">     Inderjit Dhillon, University of Texas, Austin, USA
</span><span class="comment">*</span><span class="comment">     Osni Marques, LBNL/NERSC, USA
</span><span class="comment">*</span><span class="comment">     Christof Voemel, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      INTEGER            MAXITR
      PARAMETER          ( MAXITR = 10 )
      DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, HALF
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0,
     $                     TWO = 2.0D0, THREE = 3.0D0,
     $                     FOUR = 4.0D0, HALF = 0.5D0)
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ
      INTEGER            DONE, I, IBEGIN, IDONE, IEND, II, IINDC1,
     $                   IINDC2, IINDR, IINDWK, IINFO, IM, IN, INDEIG,
     $                   INDLD, INDLLD, INDWRK, ISUPMN, ISUPMX, ITER,
     $                   ITMP1, J, JBLK, K, MINIWSIZE, MINWSIZE, NCLUS,
     $                   NDEPTH, NEGCNT, NEWCLS, NEWFST, NEWFTT, NEWLST,
     $                   NEWSIZ, OFFSET, OLDCLS, OLDFST, OLDIEN, OLDLST,
     $                   OLDNCL, P, PARITY, Q, WBEGIN, WEND, WINDEX,
     $                   WINDMN, WINDPL, ZFROM, ZTO, ZUSEDL, ZUSEDU,
     $                   ZUSEDW
      DOUBLE PRECISION   BSTRES, BSTW, EPS, FUDGE, GAP, GAPTOL, GL, GU,
     $                   LAMBDA, LEFT, LGAP, MINGMA, NRMINV, RESID,
     $                   RGAP, RIGHT, RQCORR, RQTOL, SAVGAP, SGNDEF,
     $                   SIGMA, SPDIAM, SSIGMA, TAU, TMP, TOL, ZTZ
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      DOUBLE PRECISION  <a name="DLAMCH.197"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
      DOUBLE PRECISION   DDOT, DNRM2
      EXTERNAL           DDOT, <a name="DLAMCH.199"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, DNRM2
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           DAXPY, DCOPY, <a name="DLAR1V.202"></a><a href="dlar1v.f.html#DLAR1V.1">DLAR1V</a>, <a name="DLARRB.202"></a><a href="dlarrb.f.html#DLARRB.1">DLARRB</a>, <a name="DLARRF.202"></a><a href="dlarrf.f.html#DLARRF.1">DLARRF</a>, <a name="DLASET.202"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>,
     $                   DSCAL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC ABS, DBLE, MAX, MIN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">     ..
</span>
<span class="comment">*</span><span class="comment">     The first N entries of WORK are reserved for the eigenvalues
</span>      INDLD = N+1
      INDLLD= 2*N+1
      INDWRK= 3*N+1
      MINWSIZE = 12 * N

      DO 5 I= 1,MINWSIZE
         WORK( I ) = ZERO
 5    CONTINUE

<span class="comment">*</span><span class="comment">     IWORK(IINDR+1:IINDR+N) hold the twist indices R for the
</span><span class="comment">*</span><span class="comment">     factorization used to compute the FP vector
</span>      IINDR = 0
<span class="comment">*</span><span class="comment">     IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current
</span><span class="comment">*</span><span class="comment">     layer and the one above.
</span>      IINDC1 = N
      IINDC2 = 2*N
      IINDWK = 3*N + 1

      MINIWSIZE = 7 * N
      DO 10 I= 1,MINIWSIZE
         IWORK( I ) = 0
 10   CONTINUE

      ZUSEDL = 1
      IF(DOL.GT.1) THEN
<span class="comment">*</span><span class="comment">        Set lower bound for use of Z
</span>         ZUSEDL = DOL-1
      ENDIF
      ZUSEDU = M
      IF(DOU.LT.M) THEN
<span class="comment">*</span><span class="comment">        Set lower bound for use of Z
</span>         ZUSEDU = DOU+1
      ENDIF
<span class="comment">*</span><span class="comment">     The width of the part of Z that is used
</span>      ZUSEDW = ZUSEDU - ZUSEDL + 1


      CALL <a name="DLASET.249"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, N, ZUSEDW, ZERO, ZERO,
     $                    Z(1,ZUSEDL), LDZ )

      EPS = <a name="DLAMCH.252"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Precision'</span> )
      RQTOL = TWO * EPS
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set expert flags for standard code.
</span>      TRYRQC = .TRUE.

      IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
      ELSE
<span class="comment">*</span><span class="comment">        Only selected eigenpairs are computed. Since the other evalues
</span><span class="comment">*</span><span class="comment">        are not refined by RQ iteration, bisection has to compute to full
</span><span class="comment">*</span><span class="comment">        accuracy.
</span>         RTOL1 = FOUR * EPS
         RTOL2 = FOUR * EPS
      ENDIF

<span class="comment">*</span><span class="comment">     The entries WBEGIN:WEND in W, WERR, WGAP correspond to the
</span><span class="comment">*</span><span class="comment">     desired eigenvalues. The support of the nonzero eigenvector
</span><span class="comment">*</span><span class="comment">     entries is contained in the interval IBEGIN:IEND.
</span><span class="comment">*</span><span class="comment">     Remark that if k eigenpairs are desired, then the eigenvectors
</span><span class="comment">*</span><span class="comment">     are stored in k contiguous columns of Z.
</span>
<span class="comment">*</span><span class="comment">     DONE is the number of eigenvectors already computed
</span>      DONE = 0
      IBEGIN = 1
      WBEGIN = 1
      DO 170 JBLK = 1, IBLOCK( M )
         IEND = ISPLIT( JBLK )
         SIGMA = L( IEND )
<span class="comment">*</span><span class="comment">        Find the eigenvectors of the submatrix indexed IBEGIN
</span><span class="comment">*</span><span class="comment">        through IEND.
</span>         WEND = WBEGIN - 1
 15      CONTINUE
         IF( WEND.LT.M ) THEN
            IF( IBLOCK( WEND+1 ).EQ.JBLK ) THEN
               WEND = WEND + 1
               GO TO 15
            END IF
         END IF
         IF( WEND.LT.WBEGIN ) THEN
            IBEGIN = IEND + 1
            GO TO 170
         ELSEIF( (WEND.LT.DOL).OR.(WBEGIN.GT.DOU) ) THEN
            IBEGIN = IEND + 1
            WBEGIN = WEND + 1
            GO TO 170
         END IF

<span class="comment">*</span><span class="comment">        Find local spectral diameter of the block
</span>         GL = GERS( 2*IBEGIN-1 )
         GU = GERS( 2*IBEGIN )
         DO 20 I = IBEGIN+1 , IEND
            GL = MIN( GERS( 2*I-1 ), GL )
            GU = MAX( GERS( 2*I ), GU )
 20      CONTINUE
         SPDIAM = GU - GL

<span class="comment">*</span><span class="comment">        OLDIEN is the last index of the previous block
</span>         OLDIEN = IBEGIN - 1
<span class="comment">*</span><span class="comment">        Calculate the size of the current block
</span>         IN = IEND - IBEGIN + 1
<span class="comment">*</span><span class="comment">        The number of eigenvalues in the current block
</span>         IM = WEND - WBEGIN + 1

<span class="comment">*</span><span class="comment">        This is for a 1x1 block
</span>         IF( IBEGIN.EQ.IEND ) THEN
            DONE = DONE+1
            Z( IBEGIN, WBEGIN ) = ONE
            ISUPPZ( 2*WBEGIN-1 ) = IBEGIN
            ISUPPZ( 2*WBEGIN ) = IBEGIN
            W( WBEGIN ) = W( WBEGIN ) + SIGMA
            WORK( WBEGIN ) = W( WBEGIN )
            IBEGIN = IEND + 1
            WBEGIN = WBEGIN + 1
            GO TO 170
         END IF

<span class="comment">*</span><span class="comment">        The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND)
</span><span class="comment">*</span><span class="comment">        Note that these can be approximations, in this case, the corresp.
</span><span class="comment">*</span><span class="comment">        entries of WERR give the size of the uncertainty interval.
</span><span class="comment">*</span><span class="comment">        The eigenvalue approximations will be refined when necessary as
</span><span class="comment">*</span><span class="comment">        high relative accuracy is required for the computation of the
</span><span class="comment">*</span><span class="comment">        corresponding eigenvectors.
</span>         CALL DCOPY( IM, W( WBEGIN ), 1,
     &amp;                   WORK( WBEGIN ), 1 )

<span class="comment">*</span><span class="comment">        We store in W the eigenvalue approximations w.r.t. the original
</span><span class="comment">*</span><span class="comment">        matrix T.
</span>         DO 30 I=1,IM
            W(WBEGIN+I-1) = W(WBEGIN+I-1)+SIGMA
 30      CONTINUE


<span class="comment">*</span><span class="comment">        NDEPTH is the current depth of the representation tree
</span>         NDEPTH = 0
<span class="comment">*</span><span class="comment">        PARITY is either 1 or 0
</span>         PARITY = 1
<span class="comment">*</span><span class="comment">        NCLUS is the number of clusters for the next level of the
</span><span class="comment">*</span><span class="comment">        representation tree, we start with NCLUS = 1 for the root
</span>         NCLUS = 1
         IWORK( IINDC1+1 ) = 1
         IWORK( IINDC1+2 ) = IM

<span class="comment">*</span><span class="comment">        IDONE is the number of eigenvectors already computed in the current
</span><span class="comment">*</span><span class="comment">        block
</span>         IDONE = 0
<span class="comment">*</span><span class="comment">        loop while( IDONE.LT.IM )
</span><span class="comment">*</span><span class="comment">        generate the representation tree for the current block and
</span><span class="comment">*</span><span class="comment">        compute the eigenvectors
</span>   40    CONTINUE
         IF( IDONE.LT.IM ) THEN
<span class="comment">*</span><span class="comment">           This is a crude protection against infinitely deep trees
</span>            IF( NDEPTH.GT.M ) THEN
               INFO = -2
               RETURN
            ENDIF
<span class="comment">*</span><span class="comment">           breadth first processing of the current level of the representation
</span><span class="comment">*</span><span class="comment">           tree: OLDNCL = number of clusters on current level
</span>            OLDNCL = NCLUS
<span class="comment">*</span><span class="comment">           reset NCLUS to count the number of child clusters
</span>            NCLUS = 0
<span class="comment">*</span><span class="comment">
</span>            PARITY = 1 - PARITY
            IF( PARITY.EQ.0 ) THEN
               OLDCLS = IINDC1
               NEWCLS = IINDC2
            ELSE
               OLDCLS = IINDC2
               NEWCLS = IINDC1
            END IF
<span class="comment">*</span><span class="comment">           Process the clusters on the current level
</span>            DO 150 I = 1, OLDNCL
               J = OLDCLS + 2*I
<span class="comment">*</span><span class="comment">              OLDFST, OLDLST = first, last index of current cluster.
</span><span class="comment">*</span><span class="comment">                               cluster indices start with 1 and are relative
</span><span class="comment">*</span><span class="comment">                               to WBEGIN when accessing W, WGAP, WERR, Z
</span>               OLDFST = IWORK( J-1 )
               OLDLST = IWORK( J )
               IF( NDEPTH.GT.0 ) THEN
<span class="comment">*</span><span class="comment">                 Retrieve relatively robust representation (RRR) of cluster
</span><span class="comment">*</span><span class="comment">                 that has been computed at the previous level
</span><span class="comment">*</span><span class="comment">                 The RRR is stored in Z and overwritten once the eigenvectors
</span><span class="comment">*</span><span class="comment">                 have been computed or when the cluster is refined
</span>
                  IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
<span class="comment">*</span><span class="comment">                    Get representation from location of the leftmost evalue
</span><span class="comment">*</span><span class="comment">                    of the cluster
</span>                     J = WBEGIN + OLDFST - 1
                  ELSE
                     IF(WBEGIN+OLDFST-1.LT.DOL) THEN
<span class="comment">*</span><span class="comment">                       Get representation from the left end of Z array
</span>                        J = DOL - 1
                     ELSEIF(WBEGIN+OLDFST-1.GT.DOU) THEN
<span class="comment">*</span><span class="comment">                       Get representation from the right end of Z array
</span>                        J = DOU
                     ELSE
                        J = WBEGIN + OLDFST - 1
                     ENDIF
                  ENDIF
                  CALL DCOPY( IN, Z( IBEGIN, J ), 1, D( IBEGIN ), 1 )
                  CALL DCOPY( IN-1, Z( IBEGIN, J+1 ), 1, L( IBEGIN ),
     $               1 )
                  SIGMA = Z( IEND, J+1 )

<span class="comment">*</span><span class="comment">                 Set the corresponding entries in Z to zero
</span>                  CALL <a name="DLASET.416"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, IN, 2, ZERO, ZERO,
     $                         Z( IBEGIN, J), LDZ )
               END IF

<span class="comment">*</span><span class="comment">              Compute DL and DLL of current RRR
</span>               DO 50 J = IBEGIN, IEND-1
                  TMP = D( J )*L( J )
                  WORK( INDLD-1+J ) = TMP
                  WORK( INDLLD-1+J ) = TMP*L( J )
   50          CONTINUE

               IF( NDEPTH.GT.0 ) THEN
<span class="comment">*</span><span class="comment">                 P and Q are index of the first and last eigenvalue to compute
</span><span class="comment">*</span><span class="comment">                 within the current block
</span>                  P = INDEXW( WBEGIN-1+OLDFST )
                  Q = INDEXW( WBEGIN-1+OLDLST )
<span class="comment">*</span><span class="comment">                 Offset for the arrays WORK, WGAP and WERR, i.e., th P-OFFSET
</span><span class="comment">*</span><span class="comment">                 thru' Q-OFFSET elements of these arrays are to be used.
</span>C                  OFFSET = P-OLDFST
                  OFFSET = INDEXW( WBEGIN ) - 1
<span class="comment">*</span><span class="comment">                 perform limited bisection (if necessary) to get approximate
</span><span class="comment">*</span><span class="comment">                 eigenvalues to the precision needed.
</span>                  CALL <a name="DLARRB.438"></a><a href="dlarrb.f.html#DLARRB.1">DLARRB</a>( IN, D( IBEGIN ),
     $                         WORK(INDLLD+IBEGIN-1),
     $                         P, Q, RTOL1, RTOL2, OFFSET,
     $                         WORK(WBEGIN),WGAP(WBEGIN),WERR(WBEGIN),
     $                         WORK( INDWRK ), IWORK( IINDWK ),
     $                         PIVMIN, SPDIAM, IN, IINFO )
                  IF( IINFO.NE.0 ) THEN
                     INFO = -1
                     RETURN
                  ENDIF
<span class="comment">*</span><span class="comment">                 We also recompute the extremal gaps. W holds all eigenvalues
</span><span class="comment">*</span><span class="comment">                 of the unshifted matrix and must be used for computation
</span><span class="comment">*</span><span class="comment">                 of WGAP, the entries of WORK might stem from RRRs with
</span><span class="comment">*</span><span class="comment">                 different shifts. The gaps from WBEGIN-1+OLDFST to
</span><span class="comment">*</span><span class="comment">                 WBEGIN-1+OLDLST are correctly computed in <a name="DLARRB.452"></a><a href="dlarrb.f.html#DLARRB.1">DLARRB</a>.
</span><span class="comment">*</span><span class="comment">                 However, we only allow the gaps to become greater since
</span><span class="comment">*</span><span class="comment">                 this is what should happen when we decrease WERR
</span>                  IF( OLDFST.GT.1) THEN
                     WGAP( WBEGIN+OLDFST-2 ) =
     $             MAX(WGAP(WBEGIN+OLDFST-2),
     $                 W(WBEGIN+OLDFST-1)-WERR(WBEGIN+OLDFST-1)
     $                 - W(WBEGIN+OLDFST-2)-WERR(WBEGIN+OLDFST-2) )
                  ENDIF
                  IF( WBEGIN + OLDLST -1 .LT. WEND ) THEN
                     WGAP( WBEGIN+OLDLST-1 ) =
     $               MAX(WGAP(WBEGIN+OLDLST-1),
     $                   W(WBEGIN+OLDLST)-WERR(WBEGIN+OLDLST)
     $                   - W(WBEGIN+OLDLST-1)-WERR(WBEGIN+OLDLST-1) )
                  ENDIF
<span class="comment">*</span><span class="comment">                 Each time the eigenvalues in WORK get refined, we store
</span><span class="comment">*</span><span class="comment">                 the newly found approximation with all shifts applied in W
</span>                  DO 53 J=OLDFST,OLDLST
                     W(WBEGIN+J-1) = WORK(WBEGIN+J-1)+SIGMA
 53               CONTINUE
               END IF

<span class="comment">*</span><span class="comment">              Process the current node.
</span>               NEWFST = OLDFST
               DO 140 J = OLDFST, OLDLST
                  IF( J.EQ.OLDLST ) THEN
<span class="comment">*</span><span class="comment">                    we are at the right end of the cluster, this is also the
</span><span class="comment">*</span><span class="comment">                    boundary of the child cluster
</span>                     NEWLST = J
                  ELSE IF ( WGAP( WBEGIN + J -1).GE.
     $                    MINRGP* ABS( WORK(WBEGIN + J -1) ) ) THEN
<span class="comment">*</span><span class="comment">                    the right relative gap is big enough, the child cluster
</span><span class="comment">*</span><span class="comment">                    (NEWFST,..,NEWLST) is well separated from the following
</span>                     NEWLST = J
                   ELSE
<span class="comment">*</span><span class="comment">                    inside a child cluster, the relative gap is not
</span><span class="comment">*</span><span class="comment">                    big enough.
</span>                     GOTO 140
                  END IF

<span class="comment">*</span><span class="comment">                 Compute size of child cluster found
</span>                  NEWSIZ = NEWLST - NEWFST + 1

<span class="comment">*</span><span class="comment">                 NEWFTT is the place in Z where the new RRR or the computed
</span><span class="comment">*</span><span class="comment">                 eigenvector is to be stored
</span>                  IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
<span class="comment">*</span><span class="comment">                    Store representation at location of the leftmost evalue
</span><span class="comment">*</span><span class="comment">                    of the cluster
</span>                     NEWFTT = WBEGIN + NEWFST - 1
                  ELSE
                     IF(WBEGIN+NEWFST-1.LT.DOL) THEN
<span class="comment">*</span><span class="comment">                       Store representation at the left end of Z array
</span>                        NEWFTT = DOL - 1
                     ELSEIF(WBEGIN+NEWFST-1.GT.DOU) THEN
<span class="comment">*</span><span class="comment">                       Store representation at the right end of Z array
</span>                        NEWFTT = DOU
                     ELSE
                        NEWFTT = WBEGIN + NEWFST - 1
                     ENDIF
                  ENDIF

                  IF( NEWSIZ.GT.1) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Current child is not a singleton but a cluster.
</span><span class="comment">*</span><span class="comment">                    Compute and store new representation of child.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Compute left and right cluster gap.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    LGAP and RGAP are not computed from WORK because
</span><span class="comment">*</span><span class="comment">                    the eigenvalue approximations may stem from RRRs
</span><span class="comment">*</span><span class="comment">                    different shifts. However, W hold all eigenvalues
</span><span class="comment">*</span><span class="comment">                    of the unshifted matrix. Still, the entries in WGAP
</span><span class="comment">*</span><span class="comment">                    have to be computed from WORK since the entries
</span><span class="comment">*</span><span class="comment">                    in W might be of the same order so that gaps are not
</span><span class="comment">*</span><span class="comment">                    exhibited correctly for very close eigenvalues.
</span>                     IF( NEWFST.EQ.1 ) THEN
                        LGAP = MAX( ZERO,
     $                       W(WBEGIN)-WERR(WBEGIN) - VL )
                    ELSE
                        LGAP = WGAP( WBEGIN+NEWFST-2 )
                     ENDIF
                     RGAP = WGAP( WBEGIN+NEWLST-1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Compute left- and rightmost eigenvalue of child
</span><span class="comment">*</span><span class="comment">                    to high precision in order to shift as close
</span><span class="comment">*</span><span class="comment">                    as possible and obtain as large relative gaps
</span><span class="comment">*</span><span class="comment">                    as possible
</span><span class="comment">*</span><span class="comment">
</span>                     DO 55 K =1,2
                        IF(K.EQ.1) THEN
                           P = INDEXW( WBEGIN-1+NEWFST )
                        ELSE
                           P = INDEXW( WBEGIN-1+NEWLST )
                        ENDIF
                        OFFSET = INDEXW( WBEGIN ) - 1
                        CALL <a name="DLARRB.548"></a><a href="dlarrb.f.html#DLARRB.1">DLARRB</a>( IN, D(IBEGIN),
     $                       WORK( INDLLD+IBEGIN-1 ),P,P,
     $                       RQTOL, RQTOL, OFFSET,
     $                       WORK(WBEGIN),WGAP(WBEGIN),
     $                       WERR(WBEGIN),WORK( INDWRK ),
     $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
     $                       IN, IINFO )
 55                  CONTINUE
<span class="comment">*</span><span class="comment">
</span>                     IF((WBEGIN+NEWLST-1.LT.DOL).OR.
     $                  (WBEGIN+NEWFST-1.GT.DOU)) THEN
<span class="comment">*</span><span class="comment">                       if the cluster contains no desired eigenvalues
</span><span class="comment">*</span><span class="comment">                       skip the computation of that branch of the rep. tree
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       We could skip before the refinement of the extremal
</span><span class="comment">*</span><span class="comment">                       eigenvalues of the child, but then the representation
</span><span class="comment">*</span><span class="comment">                       tree could be different from the one when nothing is
</span><span class="comment">*</span><span class="comment">                       skipped. For this reason we skip at this place.
</span>                        IDONE = IDONE + NEWLST - NEWFST + 1
                        GOTO 139
                     ENDIF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Compute RRR of child cluster.
</span><span class="comment">*</span><span class="comment">                    Note that the new RRR is stored in Z
</span><span class="comment">*</span><span class="comment">
</span>C                    <a name="DLARRF.573"></a><a href="dlarrf.f.html#DLARRF.1">DLARRF</a> needs LWORK = 2*N
                     CALL <a name="DLARRF.574"></a><a href="dlarrf.f.html#DLARRF.1">DLARRF</a>( IN, D( IBEGIN ), L( IBEGIN ),
     $                         WORK(INDLD+IBEGIN-1),
     $                         NEWFST, NEWLST, WORK(WBEGIN),
     $                         WGAP(WBEGIN), WERR(WBEGIN),
     $                         SPDIAM, LGAP, RGAP, PIVMIN, TAU,
     $                         Z(IBEGIN, NEWFTT),Z(IBEGIN, NEWFTT+1),
     $                         WORK( INDWRK ), IINFO )
                     IF( IINFO.EQ.0 ) THEN
<span class="comment">*</span><span class="comment">                       a new RRR for the cluster was found by <a name="DLARRF.582"></a><a href="dlarrf.f.html#DLARRF.1">DLARRF</a>
</span><span class="comment">*</span><span class="comment">                       update shift and store it
</span>                        SSIGMA = SIGMA + TAU
                        Z( IEND, NEWFTT+1 ) = SSIGMA
<span class="comment">*</span><span class="comment">                       WORK() are the midpoints and WERR() the semi-width
</span><span class="comment">*</span><span class="comment">                       Note that the entries in W are unchanged.
</span>                        DO 116 K = NEWFST, NEWLST
                           FUDGE =
     $                          THREE*EPS*ABS(WORK(WBEGIN+K-1))
                           WORK( WBEGIN + K - 1 ) =
     $                          WORK( WBEGIN + K - 1) - TAU
                           FUDGE = FUDGE +
     $                          FOUR*EPS*ABS(WORK(WBEGIN+K-1))
<span class="comment">*</span><span class="comment">                          Fudge errors
</span>                           WERR( WBEGIN + K - 1 ) =
     $                          WERR( WBEGIN + K - 1 ) + FUDGE
<span class="comment">*</span><span class="comment">                          Gaps are not fudged. Provided that WERR is small
</span><span class="comment">*</span><span class="comment">                          when eigenvalues are close, a zero gap indicates
</span><span class="comment">*</span><span class="comment">                          that a new representation is needed for resolving
</span><span class="comment">*</span><span class="comment">                          the cluster. A fudge could lead to a wrong decision
</span><span class="comment">*</span><span class="comment">                          of judging eigenvalues 'separated' which in
</span><span class="comment">*</span><span class="comment">                          reality are not. This could have a negative impact
</span><span class="comment">*</span><span class="comment">                          on the orthogonality of the computed eigenvectors.
</span> 116                    CONTINUE

                        NCLUS = NCLUS + 1
                        K = NEWCLS + 2*NCLUS
                        IWORK( K-1 ) = NEWFST
                        IWORK( K ) = NEWLST
                     ELSE
                        INFO = -2
                        RETURN
                     ENDIF
                  ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Compute eigenvector of singleton
</span><span class="comment">*</span><span class="comment">
</span>                     ITER = 0
<span class="comment">*</span><span class="comment">
</span>                     TOL = FOUR * LOG(DBLE(IN)) * EPS
<span class="comment">*</span><span class="comment">
</span>                     K = NEWFST
                     WINDEX = WBEGIN + K - 1
                     WINDMN = MAX(WINDEX - 1,1)
                     WINDPL = MIN(WINDEX + 1,M)
                     LAMBDA = WORK( WINDEX )
                     DONE = DONE + 1
<span class="comment">*</span><span class="comment">                    Check if eigenvector computation is to be skipped
</span>                     IF((WINDEX.LT.DOL).OR.
     $                  (WINDEX.GT.DOU)) THEN
                        ESKIP = .TRUE.
                        GOTO 125
                     ELSE
                        ESKIP = .FALSE.
                     ENDIF
                     LEFT = WORK( WINDEX ) - WERR( WINDEX )
                     RIGHT = WORK( WINDEX ) + WERR( WINDEX )
                     INDEIG = INDEXW( WINDEX )
<span class="comment">*</span><span class="comment">                    Note that since we compute the eigenpairs for a child,
</span><span class="comment">*</span><span class="comment">                    all eigenvalue approximations are w.r.t the same shift.
</span><span class="comment">*</span><span class="comment">                    In this case, the entries in WORK should be used for
</span><span class="comment">*</span><span class="comment">                    computing the gaps since they exhibit even very small
</span><span class="comment">*</span><span class="comment">                    differences in the eigenvalues, as opposed to the
</span><span class="comment">*</span><span class="comment">                    entries in W which might &quot;look&quot; the same.
</span>
                     IF( K .EQ. 1) THEN
<span class="comment">*</span><span class="comment">                       In the case RANGE='I' and with not much initial
</span><span class="comment">*</span><span class="comment">                       accuracy in LAMBDA and VL, the formula
</span><span class="comment">*</span><span class="comment">                       LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
</span><span class="comment">*</span><span class="comment">                       can lead to an overestimation of the left gap and
</span><span class="comment">*</span><span class="comment">                       thus to inadequately early RQI 'convergence'.
</span><span class="comment">*</span><span class="comment">                       Prevent this by forcing a small left gap.
</span>                        LGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
                     ELSE
                        LGAP = WGAP(WINDMN)
                     ENDIF
                     IF( K .EQ. IM) THEN
<span class="comment">*</span><span class="comment">                       In the case RANGE='I' and with not much initial
</span><span class="comment">*</span><span class="comment">                       accuracy in LAMBDA and VU, the formula
</span><span class="comment">*</span><span class="comment">                       can lead to an overestimation of the right gap and
</span><span class="comment">*</span><span class="comment">                       thus to inadequately early RQI 'convergence'.
</span><span class="comment">*</span><span class="comment">                       Prevent this by forcing a small right gap.
</span>                        RGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
                     ELSE
                        RGAP = WGAP(WINDEX)
                     ENDIF
                     GAP = MIN( LGAP, RGAP )
                     IF(( K .EQ. 1).OR.(K .EQ. IM)) THEN
<span class="comment">*</span><span class="comment">                       The eigenvector support can become wrong
</span><span class="comment">*</span><span class="comment">                       because significant entries could be cut off due to a
</span><span class="comment">*</span><span class="comment">                       large GAPTOL parameter in LAR1V. Prevent this.
</span>                        GAPTOL = ZERO
                     ELSE
                        GAPTOL = GAP * EPS
                     ENDIF
                     ISUPMN = IN
                     ISUPMX = 1
<span class="comment">*</span><span class="comment">                    Update WGAP so that it holds the minimum gap
</span><span class="comment">*</span><span class="comment">                    to the left or the right. This is crucial in the
</span><span class="comment">*</span><span class="comment">                    case where bisection is used to ensure that the
</span><span class="comment">*</span><span class="comment">                    eigenvalue is refined up to the required precision.
</span><span class="comment">*</span><span class="comment">                    The correct value is restored afterwards.
</span>                     SAVGAP = WGAP(WINDEX)
                     WGAP(WINDEX) = GAP
<span class="comment">*</span><span class="comment">                    We want to use the Rayleigh Quotient Correction
</span><span class="comment">*</span><span class="comment">                    as often as possible since it converges quadratically
</span><span class="comment">*</span><span class="comment">                    when we are close enough to the desired eigenvalue.
</span><span class="comment">*</span><span class="comment">                    However, the Rayleigh Quotient can have the wrong sign
</span><span class="comment">*</span><span class="comment">                    and lead us away from the desired eigenvalue. In this
</span><span class="comment">*</span><span class="comment">                    case, the best we can do is to use bisection.
</span>                     USEDBS = .FALSE.
                     USEDRQ = .FALSE.
<span class="comment">*</span><span class="comment">                    Bisection is initially turned off unless it is forced
</span>                     NEEDBS =  .NOT.TRYRQC
 120                 CONTINUE
<span class="comment">*</span><span class="comment">                    Check if bisection should be used to refine eigenvalue
</span>                     IF(NEEDBS) THEN
<span class="comment">*</span><span class="comment">                       Take the bisection as new iterate
</span>                        USEDBS = .TRUE.
                        ITMP1 = IWORK( IINDR+WINDEX )
                        OFFSET = INDEXW( WBEGIN ) - 1
                        CALL <a name="DLARRB.703"></a><a href="dlarrb.f.html#DLARRB.1">DLARRB</a>( IN, D(IBEGIN),
     $                       WORK(INDLLD+IBEGIN-1),INDEIG,INDEIG,
     $                       ZERO, TWO*EPS, OFFSET,
     $                       WORK(WBEGIN),WGAP(WBEGIN),
     $                       WERR(WBEGIN),WORK( INDWRK ),
     $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
     $                       ITMP1, IINFO )
                        IF( IINFO.NE.0 ) THEN
                           INFO = -3
                           RETURN
                        ENDIF
                        LAMBDA = WORK( WINDEX )
<span class="comment">*</span><span class="comment">                       Reset twist index from inaccurate LAMBDA to
</span><span class="comment">*</span><span class="comment">                       force computation of true MINGMA
</span>                        IWORK( IINDR+WINDEX ) = 0
                     ENDIF
<span class="comment">*</span><span class="comment">                    Given LAMBDA, compute the eigenvector.
</span>                     CALL <a name="DLAR1V.720"></a><a href="dlar1v.f.html#DLAR1V.1">DLAR1V</a>( IN, 1, IN, LAMBDA, D( IBEGIN ),
     $                    L( IBEGIN ), WORK(INDLD+IBEGIN-1),
     $                    WORK(INDLLD+IBEGIN-1),
     $                    PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
     $                    .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
     $                    IWORK( IINDR+WINDEX ), ISUPPZ( 2*WINDEX-1 ),
     $                    NRMINV, RESID, RQCORR, WORK( INDWRK ) )
                     IF(ITER .EQ. 0) THEN
                        BSTRES = RESID
                        BSTW = LAMBDA
                     ELSEIF(RESID.LT.BSTRES) THEN
                        BSTRES = RESID
                        BSTW = LAMBDA
                     ENDIF
                     ISUPMN = MIN(ISUPMN,ISUPPZ( 2*WINDEX-1 ))
                     ISUPMX = MAX(ISUPMX,ISUPPZ( 2*WINDEX ))
                     ITER = ITER + 1

<span class="comment">*</span><span class="comment">                    sin alpha &lt;= |resid|/gap
</span><span class="comment">*</span><span class="comment">                    Note that both the residual and the gap are
</span><span class="comment">*</span><span class="comment">                    proportional to the matrix, so ||T|| doesn't play
</span><span class="comment">*</span><span class="comment">                    a role in the quotient
</span>
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Convergence test for Rayleigh-Quotient iteration
</span><span class="comment">*</span><span class="comment">                    (omitted when Bisection has been used)
</span><span class="comment">*</span><span class="comment">
</span>                     IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT.
     $                    RQTOL*ABS( LAMBDA ) .AND. .NOT. USEDBS)
     $                    THEN
<span class="comment">*</span><span class="comment">                       We need to check that the RQCORR update doesn't
</span><span class="comment">*</span><span class="comment">                       move the eigenvalue away from the desired one and
</span><span class="comment">*</span><span class="comment">                       towards a neighbor. -&gt; protection with bisection
</span>                        IF(INDEIG.LE.NEGCNT) THEN
<span class="comment">*</span><span class="comment">                          The wanted eigenvalue lies to the left
</span>                           SGNDEF = -ONE
                        ELSE
<span class="comment">*</span><span class="comment">                          The wanted eigenvalue lies to the right
</span>                           SGNDEF = ONE
                        ENDIF
<span class="comment">*</span><span class="comment">                       We only use the RQCORR if it improves the
</span><span class="comment">*</span><span class="comment">                       the iterate reasonably.
</span>                        IF( ( RQCORR*SGNDEF.GE.ZERO )
     $                       .AND.( LAMBDA + RQCORR.LE. RIGHT)
     $                       .AND.( LAMBDA + RQCORR.GE. LEFT)
     $                       ) THEN
                           USEDRQ = .TRUE.
<span class="comment">*</span><span class="comment">                          Store new midpoint of bisection interval in WORK
</span>                           IF(SGNDEF.EQ.ONE) THEN
<span class="comment">*</span><span class="comment">                             The current LAMBDA is on the left of the true
</span><span class="comment">*</span><span class="comment">                             eigenvalue
</span>                              LEFT = LAMBDA
<span class="comment">*</span><span class="comment">                             We prefer to assume that the error estimate
</span><span class="comment">*</span><span class="comment">                             is correct. We could make the interval not
</span><span class="comment">*</span><span class="comment">                             as a bracket but to be modified if the RQCORR
</span><span class="comment">*</span><span class="comment">                             chooses to. In this case, the RIGHT side should
</span><span class="comment">*</span><span class="comment">                             be modified as follows:
</span><span class="comment">*</span><span class="comment">                              RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
</span>                           ELSE
<span class="comment">*</span><span class="comment">                             The current LAMBDA is on the right of the true
</span><span class="comment">*</span><span class="comment">                             eigenvalue
</span>                              RIGHT = LAMBDA
<span class="comment">*</span><span class="comment">                             See comment about assuming the error estimate is
</span><span class="comment">*</span><span class="comment">                             correct above.
</span><span class="comment">*</span><span class="comment">                              LEFT = MIN(LEFT, LAMBDA + RQCORR)
</span>                           ENDIF
                           WORK( WINDEX ) =
     $                       HALF * (RIGHT + LEFT)
<span class="comment">*</span><span class="comment">                          Take RQCORR since it has the correct sign and
</span><span class="comment">*</span><span class="comment">                          improves the iterate reasonably
</span>                           LAMBDA = LAMBDA + RQCORR
<span class="comment">*</span><span class="comment">                          Update width of error interval
</span>                           WERR( WINDEX ) =
     $                             HALF * (RIGHT-LEFT)
                        ELSE
                           NEEDBS = .TRUE.
                        ENDIF
                        IF(RIGHT-LEFT.LT.RQTOL*ABS(LAMBDA)) THEN
<span class="comment">*</span><span class="comment">                             The eigenvalue is computed to bisection accuracy
</span><span class="comment">*</span><span class="comment">                             compute eigenvector and stop
</span>                           USEDBS = .TRUE.
                           GOTO 120
                        ELSEIF( ITER.LT.MAXITR ) THEN
                           GOTO 120
                        ELSEIF( ITER.EQ.MAXITR ) THEN
                           NEEDBS = .TRUE.
                           GOTO 120
                        ELSE
                           INFO = 5
                           RETURN
                        END IF
                     ELSE
                        STP2II = .FALSE.
        IF(USEDRQ .AND. USEDBS .AND.
     $                     BSTRES.LE.RESID) THEN
                           LAMBDA = BSTW
                           STP2II = .TRUE.
                        ENDIF
                        IF (STP2II) THEN
<span class="comment">*</span><span class="comment">                          improve error angle by second step
</span>                           CALL <a name="DLAR1V.820"></a><a href="dlar1v.f.html#DLAR1V.1">DLAR1V</a>( IN, 1, IN, LAMBDA,
     $                          D( IBEGIN ), L( IBEGIN ),
     $                          WORK(INDLD+IBEGIN-1),
     $                          WORK(INDLLD+IBEGIN-1),
     $                          PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
     $                          .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
     $                          IWORK( IINDR+WINDEX ),
     $                          ISUPPZ( 2*WINDEX-1 ),
     $                          NRMINV, RESID, RQCORR, WORK( INDWRK ) )
                        ENDIF
                        WORK( WINDEX ) = LAMBDA
                     END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Compute FP-vector support w.r.t. whole matrix
</span><span class="comment">*</span><span class="comment">
</span>                     ISUPPZ( 2*WINDEX-1 ) = ISUPPZ( 2*WINDEX-1 )+OLDIEN
                     ISUPPZ( 2*WINDEX ) = ISUPPZ( 2*WINDEX )+OLDIEN
                     ZFROM = ISUPPZ( 2*WINDEX-1 )
                     ZTO = ISUPPZ( 2*WINDEX )
                     ISUPMN = ISUPMN + OLDIEN
                     ISUPMX = ISUPMX + OLDIEN
<span class="comment">*</span><span class="comment">                    Ensure vector is ok if support in the RQI has changed
</span>                     IF(ISUPMN.LT.ZFROM) THEN
                        DO 122 II = ISUPMN,ZFROM-1
                           Z( II, WINDEX ) = ZERO
 122                    CONTINUE
                     ENDIF
                     IF(ISUPMX.GT.ZTO) THEN
                        DO 123 II = ZTO+1,ISUPMX
                           Z( II, WINDEX ) = ZERO
 123                    CONTINUE
                     ENDIF
                     CALL DSCAL( ZTO-ZFROM+1, NRMINV,
     $                       Z( ZFROM, WINDEX ), 1 )
 125                 CONTINUE
<span class="comment">*</span><span class="comment">                    Update W
</span>                     W( WINDEX ) = LAMBDA+SIGMA
<span class="comment">*</span><span class="comment">                    Recompute the gaps on the left and right
</span><span class="comment">*</span><span class="comment">                    But only allow them to become larger and not
</span><span class="comment">*</span><span class="comment">                    smaller (which can only happen through &quot;bad&quot;
</span><span class="comment">*</span><span class="comment">                    cancellation and doesn't reflect the theory
</span><span class="comment">*</span><span class="comment">                    where the initial gaps are underestimated due
</span><span class="comment">*</span><span class="comment">                    to WERR being too crude.)
</span>                     IF(.NOT.ESKIP) THEN
                        IF( K.GT.1) THEN
                           WGAP( WINDMN ) = MAX( WGAP(WINDMN),
     $                          W(WINDEX)-WERR(WINDEX)
     $                          - W(WINDMN)-WERR(WINDMN) )
                        ENDIF
                        IF( WINDEX.LT.WEND ) THEN
                           WGAP( WINDEX ) = MAX( SAVGAP,
     $                          W( WINDPL )-WERR( WINDPL )
     $                          - W( WINDEX )-WERR( WINDEX) )
                        ENDIF
                     ENDIF
                     IDONE = IDONE + 1
                  ENDIF
<span class="comment">*</span><span class="comment">                 here ends the code for the current child
</span><span class="comment">*</span><span class="comment">
</span> 139              CONTINUE
<span class="comment">*</span><span class="comment">                 Proceed to any remaining child nodes
</span>                  NEWFST = J + 1
 140           CONTINUE
 150        CONTINUE
            NDEPTH = NDEPTH + 1
            GO TO 40
         END IF
         IBEGIN = IEND + 1
         WBEGIN = WEND + 1
 170  CONTINUE
<span class="comment">*</span><span class="comment">
</span>
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DLARRV.894"></a><a href="dlarrv.f.html#DLARRV.1">DLARRV</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
